Tumor growth instability and the onset of invasion 
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Motivated by experimental observations, we develop a mathematical model of chemotactically directed tumor 
growth. We present an analytical study of the model as well as a numerical one. The mathematical analysis 
shows that: (i) tumor cell proliferation by itself cannot generate the invasive branching behaviour observed 
experimentally, (ii) heterotype chemotaxis provides an instability mechanism that leads to the onset of tumor 
invasion and (iii) homotype chemotaxis does not provide such an instability mechanism but enhances the mean 
speed of the tumor surface. The numerical results not only support the assumptions needed to perform the 
mathematical analysis but they also provide evidence of (i), (ii) and (iii). Finally, both the analytical study and 
the numerical work agree with the experimental phenomena. 
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I. INTRODUCTION 

Experiments have shown that a variety of tumor cells pro- 
duce both protein growth factors and their corresponding re- 
ceptors, enabling a mechanism termed autocrine or paracrine, 
if the stimulus not only affects the source but also its bystander 
cells Such polypeptide growth factors can, for example, 
stimulate tumor cell growth and invasion, such as in the case 
of hepatocyte growth factors 01 and epidermal growth fac- 
tor 1 3], or induce tumor angiogenesis (through secretion of 
vascular endothelial growth factor, VEGF B). The biologi- 
cal evidence supporting these paracrine/autocrine loops sug- 
gests that such signaling factors have significance for cell- 
cell interaction |5]. Depending on the cancer type, its char- 
acteristic features include a combination of rapid volumetric 
growth and genetic/epigenetic heterogeneity, as well as ex- 
tensive tissue invasion with both local and distant dissemi- 
nation. Such tumor cell motility has been intensely investi- 
gated and found to be guided by diffusive chemical gradients, 
a process called chemotaxis, e.g., |6]. Since cell signaling 
and information processing on the microscopic scale should 
also determine the emergence of both multicellular patterns 
and macroscopic disease dynamics, it is intriguing to charac- 
terize the relationship between environmental stimuli and the 
cell-signaling code they trigger. 

In a recent paper, Sander and Deisboeck |7] showed that a 
combined heterotype and homotype diffusive chemical signal 
can yield invasive cell branching patterns seen in microscopic 
brain tumor experiments H by means of a discrete model, for 
some specific form of the interactions. There is a long history 
on the study of this type of assay ||4j,|9||. There are also quite 
a few mathematical models that address the issue of tumor 
growth and cell migration O EH E3 [3 EE C3]. 
In Ref. 0] the authors carry out a linear stability analy- 
sis from the steady state and show that both homo and het- 
erotype chemotaxis are required for the development of inva- 
sive branching behaviour. Employing an improved version of 



our previously developed reaction-diffusion model fiUl . we 
now specifically investigate the relationship between an ex- 
trinsic nutrient signal, heterotype chemoattractant Q, and the 
homotype soluble signal C, produced by the tumor cells them- 
selves. The underlying oncology concept is that in the process 
of spatio-temporal tumor expansion, C functions as a guiding 
cue for mobile cancer cells, directing the trailing ones towards 
sites of higher Q concentration and thus avoiding tissue areas 
with low or decreasing density of Q. In a sense, the dynami- 
cally changing C profile, secreted by the tumor cells, encodes 
the underlying Q map, which in itself represents a particular 
tissue environment. However, this picture is difficult to quan- 
titatively assess with conventional experimental assays, and it 
has not yet been theoretically demonstrated in a clear-cut way. 

In this paper, we are therefore particularly interested in how 
these mechanisms, homo and heterotype chemotaxis, compete 
and/or cooperate in the formation of tumor branching struc- 
tures and how the mathematical reaction terms must be cho- 
sen in order to reproduce, with some degree of universality, 
the experimentally observed patterns. Thus, we will study 
the impact of a fixed extrinsic nutrient source, and how the 
interplay among nutritive, mechanical and chemical proper- 
ties support the principle of least resistance, most attraction 
for spatio-temporal tumor cell expansion |8]. The analytical 
and numerical results provide insight as to the cells' ability 
to readily modulate C, as well as, more generally, to the im- 
portance of paracrine growth factors for information transfer 
in multicellular biosystems. We have kept the mathematical 
model as general as possible in order to understand the essen- 
tial features of the time evolution of the system. 

The report of our results is organized as follows. We de- 
scribe our reaction-diffusion mathematical model in Sec. 
where a detailed discussion of each of the involved mecha- 
nisms is given separately. Section |in]reports general analyti- 
cal results for the model of Sec.|n] We numerically check the 
validity of our analytical results in Sec. II VI Finally, we con- 
clude in Sec. |V]with a discussion of our results and provide 
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a picture of chemotactic cell invasion as triggered by tumor 
growth instability. 

n. MATHEMATICAL MODEL 

Tumor expansion is a multi-step process that involves, in a 
non-trivial fashion, several mechanisms of progression. Here, 
we concentrate on tumor cell proliferation and chemotacti- 
cally guided invasion, induced by both a diffusive heterotype 
and a homo type attractant (produced by the very tumor cells). 

A. Extracellular matrix gel (matrigel) 

The experimental setting, which we have modeled here, 
consisted of a multicellular tumor spheroid (MTS) embed- 
ded in a tissue culture medium-enriched extracellular matrix 
(ECM) gel, Matrigel [jj. We consider M = M(x, t), the av- 
erage density field of the gel matrix as a function of space x 
and time t. The role of M is twofold. From a nutrient perspec- 
tive, the tumor cells (whose concentration will be described by 
the density field U(x.,t) hereafter) metabolize M and hence 
they are able to proliferate. Besides, from a mechanical per- 
spective, the solid gel matrix has an impact on tumor cell mo- 
bility, i.e., it confines the tumor cells and so they are guided 
by least or lesser resistance areas throughout the ECM here 
in vitro, or, in vivo, by the distinct mechanical properties of 
the surrounding tissue. Therefore, at a particular site, more 
M can sustain a higher concentration of tumor cells, which in 
turn will metabolize more of the nourishing gel medium and 
thus, over time, will lead to an on-site reduction of the ECM 
matrix' mechanical resistance, which had initially hindered 
cell motility. 

Mathematically, we assume that the consumption rate of M 
by U, which we denote by Rm{M, U), grows monotonically 
with both variables, and is a non-negative function. Later, we 
will make some assumptions regarding the mechanical impact 
of M on the diffusivity of tumor cells and of both heterotype 
and homotype chemoattractants. 

For the sake of generality, we assume that the matrigel 
medium can diffuse, with constant diffusivity The order 
of magnitude of [Im depends on the specific type of medium 
under consideration (see Sec. IIII Al for further discussion on 
this issue). In summary, we can write the following equation 
for the matrigel M 

d t M = /i M V 2 M - X m Rm(M, U), (1) 

where Am is the inverse of the characteristic time of the M 
consumption process. Equation Q reflects the fact that the 
matrigel nutrient is metabolized by the tumor cells and not 
replenished. 

B. Heterotype chemoattractant 

Chemotaxis can be generally defined as motility induced 
and guided by a concentration gradient. As in our previous 



model Hill , the heterotype chemoattractant represents nutri- 
ents diffusing from a source, e.g., in vivo, a blood vessel, and 
as such is what should guide both on-site cell proliferation and 
the onset of invasion. Chemotaxis has been extensively stud- 
ied in the literature 1 19j|20J]. It is generically assumed that the 
chemotactic flux takes the form 

(Q,M)UVQ. 

Note that this flux is proportional to the tumor cell concentra- 
tion U. The function Xq(Q> M) is usually called the chemo- 
tactic sensitivity, and is a positive decreasing (or at least con- 
stant) function of both arguments Q and M . The explicit de- 
pendence of xq on M reflects the effect of the mechanical 
pressure of the underlying medium, that constrains both tu- 
mor cell and heterotype chemoattractant movement. Besides, 
tumor cells digest chemoattractant molecules, so as the former 
move towards a positive gradient of Q, the concentration of Q 
diminishes. This reduction is governed by the reaction term 
R Q (Q,U). 

Combining these ideas, we obtain the following equation 
for the heterotype chemoattractant field density Q: 

dtQ = V(mq(M)VQ) - a Q R Q {Q, U), (2) 

where ciq is the inverse of the characteristic time of the Q 
consumption process. 

In the experimental in vitro setting modeled here, the het- 
erotype chemoattractant was supplied externally. Acknowl- 
edging that the original experimental setting |8] used a non- 
replenished nutrient source, here, for simplicity, we model a 
replenished source of Q and equation (|2j thus has to be sup- 
plemented accordingly. This can be modeled by means of the 
following boundary condition: 

Q(x = L,t) = Q , (3) 

for one-dimensional systems, where L is the size of the sys- 
tem, and 

Q(x = L x ,y,t)=Q , (4) 

for two-dimensional ones, where L x is the horizontal size of 
the system and L y the vertical one. 

C. Homotype chemoattractant 

Tumor cells have been shown to produce protein growth 
factors such as the transforming growth factor alpha (TGF- 
a). These growth factors can affect the tumor producing cell 
itself, hence generate an "autocrine" feedback loop, as well 
as bystander cells, an effect called "paracrine" 12111 . In the 
following, we refer to this soluble chemical effector, as ho- 
motype chemoattractant and we denote by C its density field. 

Since an ever growing population of tumor cells digests 
more Q, the homotype chemoattractant C may take over at 
some point as guidance cue in the regions with low Q concen- 
tration. First, we assume that the homotype chemoattractant 
is both released and internalized, or (for the purposes here), 
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taken up or consumed by the tumor cells. The latter is based 
on a ligand-receptor interaction and thus on internalization of 
the class of protein growth factors, which C represents. Note, 
that if all cells produce C, a cell close to the main tumor would 
be less inclined to move away from it, since it is close to a 
large basin of C. One could tune the production rate of C in 
such a way as to ensure that only the density profile of C near 
the tumor surface has an impact on the "decision" of a tumor 
cell to stay (proliferate) or to start moving (invasion). This 
effect should have an impact on the tumor cell density U. 

The above discussion implies that C is also chemotactic for 
U. The main difference with the heterotype chemoattractant 
Q is that the homotype chemoattractant C is produced (and 
consumed) by the tumor cells. We denote by R@ (M, U) and 

Rq\c, U) the production and di gest ion rates of C, respec- 
tively. Earlier studies have shown l'22L l23l 12411 that eventually 
a central dead area develops due to the lack of nutrients in- 
side the tumor spheroid, which in turn leads to the release of 
growth inhibitory factors from the dying cells |25[. A full 
consideration of the development of such a necrotic core is 
out of the scope of this paper |26]. We consider the existence 
of this "dead area" in the dependence of R p c on the matrigel 
density M. This dependence means that the production of C 
is enhanced where M is high (outside the main tumor, as in- 
side the tumor the matrigel has been degraded) and therefore, 
the production of C is maximized for reactive tumor cells, 
i.e., surface tumor cells outside the necrotic core of the tu- 
mor |8]. We also assume that xq{Q,M) is typically larger 
than the chemotactic sensitivity of C, xc{C 1 M), for the con- 
centrations involved in the problem and that xc{C,M) de- 
pends both on C and M. The explicit dependence of xc 
on M reflects the fact that the matrigel constrains homotype 
chemoattractant movement as well. Finally, C also diffuses 
with diffusion coefficient pc(M). 

In summary, the homotype chemoattractant obeys the fol- 
lowing equation: 

8 t C = V(mc(M)VC) + a c R ( g\M, U) ~ a c R [ c\c, U), 

(5) 

where ac and ac are the inverse of the characteristic time of 
the C production and consumption process, respectively. 

D. Tumor cells 

The global nutrient density available to tumor cells is pro- 
portional to the medium density M. We then assume that 
tumor cells proliferate with a rate that is proportional to the 
rate of consumption of M. Moreover, we consider that tumor 
cells diffuse with a diffusion constant pu, and that nu de- 
pends on M to reflect the mechanical pressure of the matrigel 
M. As tumor cells are much larger than the chemoattractant 
molecules we have pQ > pc > nu- As we discussed above, 
tumor cells move towards positive gradients of both hetero 
and homotype chemoattractants, so we can write 

d t U = V(pu(M)VU) -V(U XQ {Q,M)VQ) 

-V(U X c(C, M)VC) + \uRm(M, U), (6) 



where Xu is the inverse of the characteristic time of tumor pro- 
liferation. Equations Q-ljfJ constitute our reaction-diffusion 
tumor growth mathematical model. 

III. ANALYTICAL STUDY 

As we have stated above, the precise relevance of each of 
the factors summarized in Sec. [I]] is not clearly understood. 
Partly, this is due to the complex interaction amongst these 
factors but, mainly, due to the lack of a systematic analytical 
and numerical analysis of each individual mechanism operat- 
ing in the full system. 

In this section, we analyze in detail every mechanism in- 
volved in tumor growth to determine which conditions trig- 
ger the formation of invasive branches, and how the inter- 
play among those mechanisms allows these branches to be 
sustained in time. 



A. Growth due to cell proliferation 

Consider the subsystem of equations formed by Eq. Q and 
Eq. with xq = Xc — 0, defined in a d-dimensional vol- 
ume V. The tumor and nutrient particles are confined into the 
system and so the flux of material through the boundaries van- 
ishes. This physical constraint introduces a conservation law 
in the problem, namely, 

- J dx(X M U + X a M) = 0, (7) 

independently of the precise functional form of the reaction 
term R M {M, U). 

In the absence of diffusion (pm = Hii = 0) this conserva- 
tion law provides the following closed relation 

M = K — \m/XuU, (8) 

where K depends on the initial conditions of M and U . If the 
diffusion coefficients do not vanish, we cannot obtain such 
closed relation between U and M , except in some simpler 
cases (related to the geometry of the volume and the initial 
conditions). With no loss of generality, we restrict ourselves to 
one and two-dimensional systems, V — L and V = L x x L y , 
respectively, and initial conditions such that M (x, t = Q) = 
where U(x, t = 0) ^ and U(x, t = 0) = where M (x, t = 
0) 7^ 0. Then, we can obtain a relation similar to Eq. ©. 
Physically, this means that initially the nutrient surrounds the 
implanted tumor. In this case, Eq. (|8jl is valid after a small 
transient time (see Sec. II Vt . although the shape of the front 
(i.e., surface of the tumor) changes slightly. However, as we 
are interested in the case where tumor cells diffuse slowly, this 
front will be assumed to be sharp, and hence its exact shape 
is not relevant for our discussion below fl3ll . Thereby, we 
simply get 

dtU = WipuiK-XM/XuU^^+XuRuiK-XM/XuU, U). 

(9) 
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Tumor cells digest the matrigel nutrient when they are in di- 
rect contact. Thus, the reaction term Rm (M, U) must van- 
ish when any of its arguments does. The simplest choice of 
such a reaction term is Rm(M, U) = MU . Although other 
choices are possible (leading qualitatively to the same results), 
we consider this choice to illustrate the main properties of the 
reduced system given by Eq. l|9}. 

At this stage of the formal presentation, we need to consider 
separately the cases hm -C Hu an d Mm ^> fMj- Note that 
is M-dependent so these inequalities need to be understood in 
an average sense. 



1. hm < 

We can define a small parameter e 2 — Mat/ Hu- As the dif- 
fusion coefficient of M is so small, the evolution of the density 
field M is slow in time, and therefore, random fluctuations in 
its initial condition remain at late times. Hence, the under- 
lying medium is quenched from the perspective of the tumor 
cells. Moreover, these fluctuations take place on fast length 
scales, namely, we can write fiu(M) = /if/(x/e) [27]. There 
are many studies devoted to the propagation of fronts in het- 
erogeneous media (as is the case here for late times from the 
point of view of the tumor cells) |28|. Thus, it can be shown 
that, to leading order in e, Eq. l|9} can be assumed homoge- 
neous. Namely, we can make the substitution 

fiuWe) —> ( t— ) = Mf/ = constant + 0(e), (10) 

\/W 

where (...) denotes the average over a region of length I much 
greater than the characteristic length scale of the quenched 
fluctuations of M. This means that to lowest order we can as- 
sume a constant diffusion coefficient for U and that Eq. (|9jl be- 
comes the well-known Fisher equation |29]. Fisher's equation 
admits planar traveling front solutions, with minimal wave 
speed v Q given by 13011 

«o = 2(ii u X u K) 1 / 2 . (11) 

Moreover, any deviation from the planar front (or circular for 
two-dimensional tumors) damps out, so cell proliferation by 
itself cannot generate the branches observed in our experi- 
ments (see Fig.Q. Equation provides the mean velocity 
of the tumor whenever proliferation is the only mechanism of 
tumor growth. However, chemotaxis drives tumor cells faster 
than proliferation itself, so vq is a small quantity. Thus, we 
can infer that cell proliferation is a long time process and con- 
sequently Aa/ ~ and Xjj « 0. We, therefore, assume that 
tumor cell proliferation is much slower than the chemotacti- 
cally induced tumor cell growth (see Sec. lIII BV 

Equation dlOt is only valid to lowest order in e. Corrections 
to the leading behaviour of \iu provide also corrections to the 
velocity vq. It can be shown that [28] 

v =2( t i u X u K)V 2 (l + Z 1/2 ), (12) 

where £ is obtained from the expansion of nu(M) to first or- 
der in e, and can be understood as a quenched noise term, i.e., 



FIG 1: Depicted is an overlaid image of human U87 brain tumor 
cells which were stably transfected with a Green Fluorescent Protein 
(EGFP) Histone 2B marker for nuclei. The image is taken from a 
central cross-section of the MTS cultured in a three-dimensional ex- 
tracellular matrix (Matrigel, Becton Dickinson, USA) environment 
in vitro. Note the chain-like invasive patterns. The image was taken 
one day post transferring the MTS from liquid medium to Matrigel 
(scale bar = 100 um). 



a time independent random function 12811 . Curvature correc- 
tions to Eq. J12I give the so-called quenched Kardar-Parisi- 
Zhang equation 1 3 1 ] . Thus, the heterogeneity of the matrigel 
medium will produce rough tumor interfaces. As we will see 
below, some tumor fluctuations (large length scales) are am- 
plified by chemotaxis, so they act as initial seeds for invasive 
branches. 



2. flM > H-v 

Despite the fact that the branching morphology in Fig. ^ 
has been obtained in the case um <C Mt/> for completeness, 
we include in this section the opposite limit as well. 

In this limit the diffusion coefficient of M is so large that 
any fluctuation of M is rapidly damped out. In this case, we 
cannot clearly separate the regions where M takes its limiting 
values (0 and Mo), as was possible to do in the limit \im <C 
Hu (see Fig. |2). In this case we have to deal with the full 
nonlinear equation that includes the dependence of on M. 
Hence, the specific form of the diffusion coefficient /i[/(M) is 
required in order to fully understand the evolution of the sys- 
tem. Miiller and van Saarloos 1 32] have studied the specific 
case in which (in our notation) [iu{K — Xm/XuU) ~ U k , 
with k > 0. In such case, the gradient of the tumor field den- 
sity U at the boundary of the tumor is discontinuous. This 
could be checked experimentally in order to determine the 
qualitative form of the diffusion coefficient nu(M). 
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INSIDE 
THE TUMOR 




OUTSIDE 
THE TUMOR 



FIG. 2: M density field for the two limiting cases: /xm "C (dot- 
ted line) and /xm S> /w (dashed line). The solid line represents the 
density field of tumor cells U. 



B. Growth due to heterotype chemoattraction 

In this paper we are interested in the case where the nutrient 
medium M diffuses slowly, and so the homogenization given 
by Eq. JlOi can be assumed for all diffusion coefficients and 
chemotactic sensitivities. We, therefore, drop any dependence 
of these quantities on M. In what follows we restrict ourselves 
to a two-dimensional study, with no loss of generality (the 
three-dimensional analysis can be carried out as well). 

As we have shown, cell proliferation cannot by itself pro- 
vide invasive behaviour. The next mechanism that we must 
include in order to understand cell invasion is heterotype 
chemoattraction, where the attractant molecules are provided 
externally to the tumor. They diffuse rapidly until they reach 
the tumor boundary and then two independent events take 
place: the heterotype chemoattractant is degraded by the tu- 
mor cells and the tumor cells are (chemotactically) drifted to 
higher heterotype chemoattractant concentration gradients. 

The consumption rate of the heterotype chemoattractant, 
Rq(Q, U), cannot be arbitrarily large as it saturates for large 
values of Q. This assumption is based on the concept that 
each tumor cell carries a finite number of Q-uptaking cell re- 
ceptors, which in turn determine the cell's maximum uptake 
rate. Besides, it is also a growing function of the tumor cell 
density, U. We do not need to specify the precise mathemati- 
cal form of Rq at this point, but taking into account the above 
assumptions, we can write, without loss of generality, 



R Q (Q,U) = U"<f(Q), 

with 7 a positive constant. 

In summary, the evolution equations in this case are 



and 



d t U = fi V W 2 U - V(U XQ (Q)VQ), 



d t Q = ^ Q V 2 Q~a Q UV{Q), 



(13) 



(14) 



(15) 



where we have assumed that tumor cell proliferation is neg- 
ligible compared to chemotaxis, so we can set \m = — 
vq = 0. In this limit the dynamics of M is uncoupled from 
that of Q and U. 

Despite the fact that Eqs. (114-i and dl5> are highly nonlinear, 
due to the functions xq(Q) an d f(Q)> we can obtain useful 
information by properly rescaling space, time and both field 



/ 

X 


= xf^ 


t' 


= ta Q , 


u 


= U/U , 


q 


= Q/Qo- 



densities Q and U. Thus, considering an initial tumor con- 
centration Uq located in a bounded region of the system and a 
replenished source of heterotype chemoattractant modeled by 
Eqs. (0}, we define: 

1/2 

(16) 

(17) 
(18) 
(19) 

Eqs. dl4> and dl5l can be written (we drop primes for clarity) 
as follows 

d t u = hu/hqV 2 u - vQ /{i Q V(ux(q)Vq), (20) 
d t q = X7 2 q-u~<U^f(Qoq)/Q Ql (21) 

where x(q) is a dimensionless version of xq(Q) an d v is de- 
fined through the relation v = Xq{Q) / x(q)- 

Note that /iq is the fastest diffusivity in the problem, so we 
can define an small parameter e = hu/pq. Moreover, we 
also assume that the cross diffusion coefficient vQ§ is smaller 
than hq 13311 . With these considerations Eq. ( 1201 becomes 



d t u 



eV 2 u 



peV{ux{q)Vq), 



(22) 



with p = vQ /{j,[j. Typically, the tumor field will be con- 
stant almost everywhere except in a narrow region (boundary 
layer 1 34]). This region defines an interface between the in- 
side and the outside of the tumor. 

In fact, if we take the limit e — > (the so-called outer 
limit 1 34]) in Eq. \22\ . we have dtu = 0, and 



inside the tumor, 
outside the tumor. 



(23) 



One can see that in the outer limit the mathematical analy- 
sis of the problem is much simpler as Eq. J2 1 i reduces to the 
equations 



dtq={ 



V 2 9 - U^f(Q q)/Q 
V 2 g 



inside the tumor, 
outside the tumor. 



(24) 



In the case e ^ (the so-called inner limit 1 34]) we need to 
proceed with care: in order to solve Eq. J20i for u and Eq. Mil 
for q, we need the behavior of q exactly at the tumor boundary 
layer (interface). With this in mind, we define a new local set 
of curvilinear coordinates: a coordinate n normal to the tumor 
interface and a coordinate s tangential to it. Elementary com- 
putations (see Ref. l35ln give us the formulas for converting 
derivatives with respect to x to derivatives with respect to the 
new coordinates (n, s): 



V 2 =d„ 



(25) 



where k is the local dimensionless curvature of the tumor in- 
terface and A s is the surface Laplacian 1 35 ] . Similarly, we 
can write 



d t = d t - vd n + s t d s , 



(26) 
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where v is the normal (dimensionless) velocity of the inter- 
face. We also need the following derivative 

V(ux(g)Vq) = d n (ux(q)d n q)+R(ux{q)d n q)+V s {ux(q)V s q) 

(27) 

Our aim now is to find solutions for q and u in the tumor 
interface. We start by rescaling the normal coordinate n, in 
terms of the fast variable 77 = n/e, and time in terms of r = 
t/e. We denote by u and q the inner fields, and expand them, 
x(q) and v in a power series of e 13611 : 



which states the fact that, to first order, the tumor cells simply 
diffuse. Note that to this order uo does not depend on s, so 
that d s uo — 0, due to the boundary conditions on u. 



2. Order e 1 



To next order in e and from Eq. d33l we find, 



u(rj,s,t) = uo + cui + ..., 

q(v, s,t) = q + eq% + 

v(r),s,t) = v + evi + . . . , 

x(q) = x(q~o) + x'(qo)eqi 



(28) 
(29) 
(30) 
(31) 



where the prime in Eq. i3l\ denotes a derivative with respect 
to q. Using Eqs. i22\ . (I25> . (I26> and ( I27> . the original system 
of equations can be written as 



d vv qi — => d v qi = constant, 

due to the boundary conditions on q. 

The next order equation for u\ is given by 



d T uo + ed T ui — 



ed T q 



evid^uo + d vv u 
ed vv iii — estdgito + ekd v uo 
pd v (uox(qo)d v qo) - epd v (uix(qo)dr,qo) 
£pdr,(u x(qo)d n qi) - epd n (u x' (q )qid v qo) 
epKu x(qo)d v qo + 0(e 2 ), 
d vv qo + td m qi + £K,d v q + 0(e 2 ). 



(32) 
(33) 



1. Order e' 



We start by solving Eqs. (I32> and J33i to order e°. As we 
stated above, tumor cell proliferation is negligible when com- 
pared to chemotaxis, which means that vq = 0. Otherwise 
we have to include the reaction term Rm in the dynamical 
equation for U and then, to order e°, Eq. d32i becomes: 

d T u + ed T ui — vod^uo + tvxd^uo + ev d v ui + d vv UQ 
+ ed m ui — es t d s uo + ekd ri v,Q 

- pd n (uox(q )d v q ) - ep9r ) (Mix(*)9 7) go) 

- epd n (u a x(qo)d n qi) - epd v (u x'(qo)qid v qo) 

- £pKUox(qo)d v qo 

0(e 2 ) 



(IQ 



-u u 



XuK 



(37) 



d T u x = v 1 dy 1 uo + d m u 1 + kd 11 uo-p{d r] uo)xi.qo)d v qi- (38) 

The Fredholm alternative (or solvability condition) for ui pro- 
vides 136] 



Vi 



where 



B 



-K + pBd v qi, 

f-oo d v(d n uo) 2 x(qo) 



JZo dr](d v u ) 2 



(39) 



(40) 



Matching the inner and outer expansions yields the normal 
velocity of the tumor in the original dimensions 1 36]: 



v n = —pu K + vBXIQ ■ n, 



(41) 



where n is a unit vector normal to the tumor interface and 
directed away from the tumor. 

Chemotactically induced tumor growth requires VQ • n to 
be a positive quantity, so that v n is positive. This can be 
achieved whenever the heterotype chemoattractant concentra- 
tion grows as we move away from the tumor surface. Before 
performing the analysis of Eqs. ( 1241 and J41i to check this 
requirement, we consider two interesting features related to 
Eq. \A\\ . First of all, in the absence of chemotaxis and prolif- 
eration, for an initially circular tumor of radius Rq, Eq. ( 14 11 
reduces to 



dR 
~dt 



Pu 
' R ' 



(42) 



(34) 



Notice that in the absence of chemotaxis Eq. J34i becomes 
Fisher's equation in the set of coordinates (77, s) and with time 
given by r. 

From Eq. ( I33l > we find 



d m qo = 0. 



(35) 



We know that qo must be bounded at 77 = ±00, corresponding 
to the inside and outside of the tumor in the scaled variable 
77. This implies d n qo = and, therefore, qo — constant. 
Substituting this into Eq. d32i . we obtain 



which gives the diffusive behavior R(t) = \J Rq — 2ujjt. 
This means that the mean velocity of the tumor due to diffu- 
sion decays as t~ 1//2 . At this stage, and for an initially circular 
tumor, one could now perform the following analysis: (i) sup- 
pose the tumour continues to grow as a two-dimensional disc 
and (ii) perturb the boundary and study the development of 
instabilities. This is out of the scope of this paper and will be 
addressed in future work. Secondly, for a general tumor front 
geometry, Eq. (14- 1 i provides a critical curvature, (reminiscent 
of the classical Greenspan model llTHo 



k c = VQ ■ n, 

Pu 



(43) 



d T U0 = <9 W M , 



(36) 



such that for k < k c the tumor interface has locally a positive 
velocity, for k > k c the tumor interface has locally a negative 
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velocity and for k — k c the tumor interface has locally van- 
ishing velocity. Note that this is precisely the reason of the 
dynamical instability that guarantees the development of in- 
vasive branches. Tumor invasion will take place on the tumor 
surface wherever the local curvature is below the critical value 
K c , given by Eq. ( |43> . Our result agrees with the specific case 
considered in Ref. Il7ll . 

As a final remark concerning Eq. J41I . note that it resembles 
the equation for the velocity of a solidification front I37 U3811 . 
In that case, the local normal velocity of the front depends 
on the local curvature, k as above, but on the value of the 
field (temperature) at the boundary, instead of the value of the 
gradient, VQ, at the tumor interface. This is to be expected as 
in the tumor picture it is the heterotype chemoattractant that is 
driving the dynamics, and in the case of a solidification front, 
the dynamics of the front is linked to the temperature field 1371 
13^1 . The branching behaviour of the tumor also resembles the 
fingering instability of the Hele-Shaw problem l39ll 

We now turn to the requirement on the value of VQ • n. In 
order to check that, indeed, VQ • n (or equivalently Vg • n) 
is a positive quantity at the tumor boundary, we need to solve 
Eq. i24l supplemented by Eq. J41I . with initial condition, 



q(x,y,t = 0) = S(x - L x ), 

and boundary conditions 

q(x = 0,y,t) = 0, 
q(x = L x ,y,t) = 1. 



(44) 



(45) 
(46) 



These equations cannot be solved in general without the pre- 
cise form of the function f(Q). Nevertheless, we may assume 
that, as tumor cells degrade the heterotype chemoattractant, 
the concentration of q inside the tumor is small enough, and 
we can approximate f(Q) by a linear function f(Qoq) — Sq. 
We now proceed by solving Eq. (124b inside and outside the 
tumor (we denote the solutions, respectively, by q~ and q + ) 
and matching both solutions at the boundary, for an initially 
flat tumor, i.e., a tumor for which the radius of curvature is 
much smaller than the system lateral size L x & We first 
notice that q + is invariant under the transformation group 
(L x — x) — > e(L x — x), t — > e 2 i and q + — > e a q + . Simi- 
larly, e 



St q is invariant under the group, x — > ex, t — > e 2 t and 
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e q . Thus, it can be straightforwardly seen that 1 4 1 



q (x,V,t) = Ce 



-Si 



x/y/t 



d£,e 



(47) 



q + (x,y,t) = 1-A / dte-t '\ (48) 

Jq 

where A and C are positive constants that can be deter- 
mined by continuity of the solutions at the tumor boundary 
x = Xo(s, t). As we are interested in the gradient of the Q 
density field, we find 



d x q (x,y,t) 
d x q + (x,y,t) 



C_ 

Vt 

A 



e «e-*V4* 



-(L,-i) 2 /4i 



(49) 
(50) 



Note that Eq. J50i states that the gradient of Q is a positive 
function, and so the tumor velocity is increased by chemotaxis 
as we had anticipated. Moreover, the larger the distance from 
x to the tumor is, the larger the value of d x q + becomes and 
the larger the value of the velocity front is (see Fig.|5|l. This is 
indeed, the reason why small fluctuations on the tumor surface 
become emerging invasive branches. We will see in Sec. IIVI 
that this chemotactically enhanced velocity is also obtained 
numerically. 




TUMOR 



FIG. 3: Sketch of the heterotype chemoattractant effect on tumor cell 
invasive drift. Dashed lines represent level sets of Q and arrows rep- 
resent the local normal velocity of the tumor cells due to chemotaxis 
(the length of the arrows is proportional to the normal velocity). 



C. Growth due to homotype chemoattr action 

Finally, we consider the role of homotype chemoattractants. 
Representing protein growth factors, these chemoattractants 
are produced and internalized ll42ll . or (for the purposes here), 
consumed by the tumor cells that move towards their positive 
gradients, in a similar fashion as in the heterotype case. Yet, 
there is no wide time scale separation between tumor and ho- 
motype chemoattractant dynamics. Thus, we cannot, in gen- 
eral, neglect tumor cell proliferation when analyzing homo- 
type chemoattractant dynamics. This means that the equations 
in this case are Eqs. Q, (|5} and with \Q — 0- Performing 
a similar analysis to that of Sec. lHIBl above we find: 



v n = 2{p, u X u K) 1/2 - [MxjK + v'B'VC ■ n. 



(51) 



Despite the fact that Eq. ( 15 H is equivalent to Eq. 14 H . the 
main differences between the evolution of Q and C are due to 
the behaviour of both density fields away from the tumor in- 
terface, namely, due to the production term Rq (that is absent 
in the dynamical equation for Q) and the different boundary 
conditions (no external source for C). 

Following the same steps as those carried out in Sec. lHIBl 
we find that the outer limit yields the following equation 



d t c - 



V 2 c + rW(0, 1) -rW{c, 1) 



id) I 



inside the tumor, 
outside the tumor, 
(52) 
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with r( p ) and scaled (dimensionless) versions ofRg> and 

Rq , respectively. Note that we have considered m = 
where u = 1. The conservation relation given by Eq. (|8jl still 
holds in this case so, without loss of generality, we reduce the 
system of equations Q, © and to Eq. (jSJi and Eq. 0, with 
Xq = and M given by Eq. (|8j. 

Just as we did in Sec. 1111 Bl we must now compute the 
sign of VC • n in order to determine whether or not homo- 
type chemotaxis increases the velocity of the tumor bound- 
ary, and if it can generate a dynamical instability leading to 
a tumor branching morphology. We cannot in general find a 
transformation group under which Eq. J52b is invariant. This 
means that we cannot study homotype chemotaxis with the 
tools of the previous section. However, we can analyze homo- 
type chemoattraction by means of its homogeneous and steady 
state solutions, based on generic assumptions regarding the re- 
action terms R ( g ] [M, U) and (C, U). 

The homogeneous, steady state solutions (nullclines 1 4-3] ) 
of equations (|9jl and @ are given by the solutions of 

R M (K -X M /XuU,U) = 0, (53) 

R^\K-X M /XuU,U) = R ( c\c,U). (54) 

These nullclines yield the fixed points U\ = and U2 — 
XtjK/Xm, and C\ and C2 given implicitly by Eq. ( I54> . We 
must distinguish between the inside and the outside of the tu- 
mor when analyzing tumor growth due to homotype chemo- 
taxis. Outside the tumor U% = and there is no production 
of C (as there are no tumor cells). This implies that the corre- 
sponding fixed point value of C is then C\ = 0. On the other 
hand, inside the tumor, Eq. d54i reflects the fact that there is 
a balance between production and consumption of C. This 
means that the value of the density field C reaches the fixed 
point C2, which is a constant (equilibrium) value. This simple 
picture holds even far enough from the tumor boundary (dif- 
fusion tends to spread these uniform concentration phases). 
Therefore, let us consider a point in the C — U plane where 
U = Up and C — (point P in Fig. [4^), namely, a point 
outside the tumor boundary, with U so small that there has 
been no previous secretion of homotype chemoattractant C. 
The concentration C eventually grows (as U increases due to 
proliferation wherever M ^ 0) because the slope given by 

dC = R™(K - Xm/XuU, U) - rW(C, U) 
dU R M (K-X M /XuU,U) 

is positive 1 44] . Then, the slope decreases |45] until it reaches 
the nullcline where the slope is (point Q in Fig.|4^). Finally, 
the curve approaches the stable point (1/2,02) with infinite 
slope (point R in Fig. |4^). Thus, C grows from the outside 
of the tumor, reaches a maximum value and then decreases to 
a constant value C2 inside the tumor. This can be schemat- 
ically seen in Fig. |4j?). In other words, the density field C 
tends to grow as U decreases namely, as we move outside of 
the tumor from inside. But, as we have shown, outside the 
tumor and far enough from it C tends to 0. This means that 
C also tends to grow as we move inside of the tumor from the 
outside. Consequently, C behaves as a pulse from the con- 
stant value C% inside the tumor to C\ — outside the tumor. 



This qualitative picture is sketched in Fig.|4] where we have 
shown that the nullcline given by Eq. j54l > can have three dif- 
ferent qualitative behaviours N±, N2 and N3, that are plotted 
as well. This pulse-like structure for the density profile of C 
agrees well with the intuitive picture provided in Sec. Ill CI It 
also explains why tumor cells are chemotactically guided by 
the gradient of C. These facts will be numerically confirmed 
in Sec. llVl below. We conclude as follows: the density profile 
of C grows at the tumor boundary, which implies that the lo- 
cal normal velocity of the tumor surface increases due to the 
presence of a positive homotype chemoattractant gradient. 



a) b) 




FIG. 4: a) Qualitative behaviour of the trajectory of the density field 
C in the phase-plane. The arrows indicate the direction of time evo- 
lution. Curves iVl-iVS are the three types of nullclines that can be 
expected for our system, b) Same trajectory of the density field C 
with respect to the spatial coordinate x. Again, arrows indicate time 
evolution. 



IV. NUMERICAL STUDY 

In Sec. HJl] we have presented a general analytical frame- 
work to study tumor growth (proliferation and chemotactic 
invasion). The previous analysis has been carried out by 
means of several assumptions and limits, {e.g., vQ -C [iq, 
Am ~ 0, Xjj w 0), that have not been fully justified. In this 
section we provide numerical simulations that check the va- 
lidity of those assumptions and limits. We, thus, numerically 
solve the main differential equations of Sec. [H] following a 
similar organization to that of Sec. [HI] It is clear that in order 
to carry out a numerical study we must specify the detailed 
mathematical form of the reaction terms in our equations Q- 
l|6}. These reaction terms are chosen following the spirit of 
the oncology concept previously reported in Refs. 1 8. I18I1 . 

A. Growth due to cell proliferation 

In this section, we check the validity of Eq. as an ap- 
proximation to Eq. @. Thus, we have integrated Eqs. Q 
and l|6} with xq — Xc = and Rm(U,M) = MU and, 
independently, Eq. l|6} with M = K — Xm/XjjU . 



1. One-dimensional results 

Initially, we place a one dimensional tumor such that 
U(x,t = 0) = 1 for < x < L/2 and U(x,t = 0) = 
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elsewhere. This implies that at the initial time M (x, t = 0) = 
1 — U(x, t = 0). We have numerically solved Eqs. Q and 
with xq = Xc = and Rm(U,M) = MU and, indepen- 
dently, Eq. (jfjl with M = K — Xm/XjjU. The lattice separa- 
tion has been taken dx — 0.5, the lattice size L x = 128, the 
diffusion coefficients fiu = 0.01 and /jm = 0.0, the reaction 
rates Am = Xu = 0.1, the time step e f = 0.005, the initial 
time to = € t and the final time f/ = 50000e t . 

As can be seen in Fig. [5] after a small transient time, the ap- 
proximation is accurate enough. As we mentioned in Sec. Mil 
the shape of the tumor front changes slightly. 




X 



dx = dy = 0.5, the lattice size L x = L y = 128, the diffu- 
sion coefficients fLrj = 0.01 and \lm = 0.0, the reaction rates 
Am = Ay = 1.5, the time step et — 0.005, the initial time 
to = et, the final time tj = 5000et, and the intermediate times 
(1000, 2000, 3000, 4000)e t . 

Figure [6] displays the time evolution of the tumor surface. 
Notice how the tumor conserves during its evolution its ini- 
tial circular shape but develops a rough interface with the ma- 
trigel. 




FIG. 6: Numerical simulation of an initially circular tumor embedded 
in a matrigel medium M with slow dynamics (/xm <C p-u). Different 
curves represent different times with the initial time to corresponding 
to the inner perfect circle. 



FIG. 5: One dimensional tumor proliferation. Solid lines represent 
the U density field for the model given by Eqs. and |6j. Dashed 
lines represent U for the solution of Fisher's equation, i.e., Eq. JSJ 
with M = K - Am f\uU. 



2. Two-dimensional results 

We now consider tumor growth due to proliferation in a 
two-dimensional setting with \Q = Xc = 0. Initially, the 
matrigel density M is a random distribution to reflect the fact 
that it is an heterogeneous medium. As we are interested in 
a slowly varying nutrient, we choose e = hm/P-u = 0.01, 
where Jiu is the maximum attainable value of fiu(M). We 
assume that this value of /j,tj corresponds to the value M = 0. 
The confinement due to the matrigel is unlimited, namely, for 
large concentrations of M, tumor cells can no longer diffuse. 
Hence, we take the following diffusion coefficient 



MM) 



1 + M/M th ' 



(56) 



where M t h is a reference threshold concentration. 

We have solved Eqs. Q and with the following initial 
conditions: at time t = 0we place a circular tumor centered at 
(L x /2, Ly/2) of radius L K /4 and surrounded by a heteroge- 
neous nutrient substrate M. Thus, M(x, y, t = 0) = inside 
the initial circular tumor and M(x, y, t = 0) = 1 + y) 
elsewhere, with £(x,y) a random Gaussian distribution with 
zero mean and variance 0.2, that encodes the initial hetero- 
geneities of the matrigel. The lattice separation has been taken 



B. Growth due to heterotype chemoattraction 

The experimental branches of Fig.^have two characteris- 
tic lengths, namely, their width and their height with respect to 
the main tumor substrate. Section lTiIBI was devoted to deter- 
mine the conditions that trigger the formation of the invasive 
branches. We know that the height of the branches depends in 
a crucial way on the mathematical form of the reaction term 
Rq(Q, U) of Eq. (|2j, namely, the height depends on the con- 
centration thresholds associated with the bio-chemical reac- 
tion of Q consumption by tumor cells 0]. Thus, we choose 



aQRqiQ^U) = a Q U 



Q 



b Q 



(57) 



where ag is the inverse of the time scale of consumption and 
6q a characteristic heterotype chemoattractant concentration 
(threshold value). Note that, in the notation of Eq. H51 we 
have set 7 = 1. 



1. Two-dimensional results 

We have integrated Eqs. 0, (|6j and (|2j with Rq given by 
Eq. ( I57l > and \c = 0. At the initial time we place a circu- 
lar tumor at (L x /2, L y /2) with radius L x /A + £, where £ 
is a Gaussian random number with zero mean and variance 
L x /20. This initial condition mimics the effect of the slowly 
varying underlying substrate (matrigel) as shown in the pre- 
vious section IIV Al The lattice separation has been taken 
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dx = dy = 0.5, the lattice size L x = L y = 128, the diffusion 
coefficients ^ v = 0.001, \i Q = 10.0 and fi M = 0.00001, the 
reaction rates Xm = \u = 0.0050, q,q = 0.75, bq — 0.5, the 
chemotactic sensitivity \Q = 2-0, the time step e t = 0.005, 
the initial time to = e t and the final time tf = 20000et. 

Figure Et) displays the numerically obtained tumor when 
the replenished source of Q is placed at x = L x (i.e., the right 
hand side of the lattice). Moreover, in Fig. [7J)) we show the 
cross-section of an invasive, chemotactically induced branch 
obtained in the same simulation. Clearly, we can distinguish 
between the main tumor spheroid and a given branch. Notice 
the good agreement with the experimental results and with the 
qualitative analysis provided in Sec. UTTI 




FIG. 7: a) Numerical simulation of tumor branching induced by a 
heterotype chemotactic source located at x = L x . The plot rep- 
resents the tumor density field U(x,y,t = 20000et). b) Cross- 
section of the tumor in panel a) for different times between 2000et 
and 20000e t . 



C. Growth due to homotype chemoattraction 



same as those chosen in the previous section ITV Bl For the 
matrigel and the homotype chemoattracctant, we have chosen 
M (x, y, t = 0) = 1 - U(x, y,t = 0) and C(x, y, t = 0) = 0, 
respectively. The lattice separation has been taken dx = dy = 
0.5, the lattice size L x — L y = 128, the diffusion coefficients 
Hu — 0.01, fxc = 1-0 and ^ = 0.0, the reaction rates 
A M = Xu = 0.1, a c = 1.75, b c = 0.1, a c = 0c = 1.0, the 
chemotactic sensitivity xc = 1-0, the time step e t = 0.005, 
the initial time to — e t and the final time tf = 50000e(. 

Fig- HI displays the time evolution of the tumor profile in 
the x direction for times 10000e t , 15000e t , 20000e t , 25000e t 
and 30000e(. The dotted line in Fig.[8]displays the C profile 
at time t = 10000et, magnified by a factor of 4. In agree- 
ment with the analysis presented in section UlI CI the numer- 
ical results show that (i) there is no emergence of chemotac- 
tically induced branches, as was the case for the heterotype 
chemoattractant and (ii) the mean speed of the tumor bound- 
ary increases due to C. That is, the tumor profile follows qual- 
itatively the behaviour anticipated in Sec. lHICl 




X 



The segregation and eventual degradation of the homo- 
type chemoattractant is limited, i.e., the rates associated with 
both processes cannot be arbitrarily large. Thus, following 
Ref. (H we choose 

a c R < g ) (U) = ac- r ^- F 7, (58) 



FIG. 8: Cross-section of a tumor for the case of coupled dynamics 
between matrigel M, tumor cells U and homotype chemoattractant 
C for different times between to = e t and tf — 30000et- Solid 
lines: U subject to homotype chemotaxis; dashed lines: U subject to 
no homotype chemotaxis. Dotted line: C (times a factor of four) for 
t = lOOOOet. 



a c R { c\C, U) = a c U——. (59) 

The constants ac and ac are the inverse of the characteris- 
tic time scales of production and degradation, respectively, 
and fie and be are characteristic saturation concentrations (for 
production and degradation, respectively). Note the Eqs. d59l 
and \51\ have the same mathematical form. 

1. Two-dimensional results 

We have numerically integrated Eqs. Q, and Q with 
Xq = 0, and and given by Eqs. d58l and ( I59> . re- 
spectively. The initial conditions for U(x,y,t = 0) are the 



V. DISCUSSION AND CONCLUSIONS 

In summary we conclude that: 

1 . The matrigel M induces tumor cell proliferation. This 
growth is an overall expansion of the initial tumor that 
follows the principle of least resistance 1 18]. Moreover, 
in the case of interest here, that of a slowly diffusing 
matrigel, the tumor boundary (or surface) becomes in- 
homogeneous due to the random nature of the slowly 
varying nutrient M. It is noteworthy that the roughness 
of the tumor surface depends also on the proliferation 
rate of the cells. 

2. The homotype chemoattractant enhances the velocity of 
the tumor cells due to the increase of C at the boundary 
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of the tumor. This effect combined with cell prolifer- 
ation (due to M) would induce the onset of invasion 
of tumor cells towards regions of lower matrigel den- 
sity. We have been able to show that the secretion and 
subsequent diffusion of C catalyzes the motion of the 
tumor cells. 

3. Finally, as the heterotype chemoattractant is introduced 
into the system at a given distance from the tumor 
(in our case at the boundary), and because it diffuses, 
an initially circular tumor develops unstable invasive 
branches that move towards the source of the heterotype 
chemoattractant. These branches develop from initial 
seeds (i.e., fluctuations on the tumor boundary), which 
are due to the mechanical confinement from M, and the 
velocity enhancement due to C and Q. 

Our results are an improvement over We have not 

limited ourselves to (i) performing a linear stability analy- 
sis from the steady state solutions as was done in 0] or to 
(ii) numerically solving a simplified version of the reaction- 
diffusion equations as carried out in 1 18]. On the other hand, 
we have analytically and numerically studied the full non- 
linear problem. The results of the work presented here allow 
us to say that proliferation is a requirement for invasion. In 
fact, there can be no (chemotactic) invasion without prolif- 
eration, as proliferation due to M provides the initial seeds 
that trigger the onset of invasion. That is, the slow diffu- 
sion of M is crucial to the development of those initial seeds 
(rough tumor interface) that become invasive branches due to 
heterotype chemotactic ally induced instability. 

Admittedly, our model still does have several shortcomings 
as it inevitably has to simplify the complex biological sce- 
nario considered here. For instance, tumor cell apoptosis and 
thus, the development of a central necrotic core is currently 
not included 1 26] . Incorporating this characteristic tumor fea- 
ture would also have implications for the simulation itself. 
Specifically, detrimental byproducts released from the dying 
virtual cells would render this area "toxic", resulting in a cen- 
tral space within the growing tumor, which is not being re- 
populated by the proliferative tumor surface. In future work, 



this tumor characteristic can be implemented e.g., by some dy- 
namic, internal boundary condition within the tumor. We have 
also failed to model the finite receptor occupancy of the tumor 
cell surface. This issue is important in order to find biologi- 
cal support for implementing a maximum threshold of (e.g., 
homotype) chemoattractant uptake rate (by each tumor cell). 
On the other hand, the minimum threshold is given by the 
maximum sensitivity of the cell surface-based receptor sys- 
tem. Nonetheless, even in its present form the model already 
proves very useful for interdisciplinary cancer research as it 
provides the following, at least in part experimentally testable 
hypotheses: (i) tumor cell proliferation by itself cannot gener- 
ate the invasive branching behaviour observed experimentally, 
yet, proliferation is a requirement for invasion (ii) heterotype 
chemotaxis provides an instability mechanism that leads to the 
onset of tumor cell invasion and (iii) homotype chemotaxis 
does not provide such a mechanism but enhances the mean 
speed of the tumor surface. 

Combined with more specific experimental data, both on 
the molecular and on the microscopic scale, this ongoing work 
may therefore reveal novel and exciting insights into the role 
of tumor cell signaling and its impact on the emergence of 
multicellular patterns. 
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